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Abstract. The Kosambi-Cartan-Chern (KCC) theory represents a powerful 
mathematical method for the analysis of dynamical systems. In this approach 
one describes the evolution of a dynamical system in geometric terms, by 
considering it as a geodesic in a Finsler space. By associating a non-linear 
connection and a Berwald type connection to the dynamical system, five geo- 
metrical invariants are obtained, with the second invariant giving the Jacobi 
stability of the system. The Jacobi (in)stability is a natural generalization 
of the (in)stability of the geodesic flow on a differentiable manifold endowed 
with a metric (Riemannian or Finslerian) to the non-metric setting. In the 
present paper we review the basic mathematical formalism of the KCC the- 
ory, and present some specific applications of this method in general relativity, 
cosmology and astrophysics. In particular we investigate the Jacobi stability 
of the general relativistic static fluid sphere with a linear barotropic equation 
of state, of the vacuum in the brane world models, of a dynamical dark energy 
model, and of the Lane-Emden equation, respectively. It is shown that the Ja- 
cobi stability analysis offers a powerful and simple method for constraining the 
physical properties of different systems, described by second order differential 
equations. 
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1. Introduction 

A second order differential equation can be investigated in geometric terms by us- 
ing the general path-space theory of Kosambi-Cartan-Chern (KCC-theory) inspired 
by the geometry of a Finsler space jT5J [HI U\ ■ The KCC theory is a differential geo- 
metric theory of the variational equations for the deviation of the whole trajectory 
to nearby ones. By associating a non-linear connection and a Berwald type connec- 
tion to the differential system, five geometrical invariants are obtained. The second 
invariant gives the Jacobi stability of the system [23l EI]- The KCC theory has 
been applied for the study of different physical, biochemical or technical systems 
(see Ell EH HJ El [27] and references therein). 

The Jacobi stability of a dynamical system can be regarded as the robustness 
of the system to small perturbations of the whole trajectory [23]. This is a very 
convenable way of regarding the resistance of limit cycles to small perturbation of 
trajectories. On the other hand, we may regard the Jacobi stability for other types 
of dynamical systems (like the ones in the present paper) as the resistance of a whole 
trajectory to the onset of chaos due to small perturbations of the whole trajectory. 
This interpretation is based on the generally accepted definition of chaos, namely 
a compact manifold M on which the geodesic trajectories deviate exponentially 
fast. This is obviously related to the curvature of the base manifold (see section 
13. 1|) . The Jacobi (in)stability is a natural generalization of the (in)stability of the 
geodesic flow on a differcntiable manifold endowed with a metric (Riemannian or 
Finslerian) to the non-metric setting. In other words, we may say that Jacobi 
unstable trajectories of a dynamical system behave chaotically in the sense that 
after a finite interval of time it would be impossible to distinguish the trajectories 
that were very near each other at an initial moment. 

It is the purpose of the present paper to review the KCC theory, to analyze 
its relations with the linear Lyapunov stability analysis of dynamical systems, and 
to present a comparative study of these methods in the fields of gravitation and 
astrophysics. The Jacobi stability analysis has been extensively applied in different 
fields of science and technology, but presently its applications to the very important 
fields of gravitation, cosmology and astrophysics have been limited to the Jacobi 
stability analysis of the Lane-Emden equation 0] , and to the analysis of the stability 
of the vacuum in the brane world models |12) . Besides reviewing these examples, 
we will present in detail several new applications in gravitation and cosmology of 
the KCC theory, including the analysis of the general relativistic static fluid sphere, 
and of the interacting dark energy models. In all these cases we will also carefully 
consider the physical implications of our stability analysis. 
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The present paper is organized as follows. In Section [2] we review the basic 
mathematical concepts related to the analysis of the linear stability of the dynam- 
ical systems. The KCC theory and the Jacobi stability analysis of the dynamical 
systems are presented in Section[3l The case of the Newtonian polytropes is consid- 
ered in Section[4] The Jacobi stability analysis of the general relativistic fluid sphere 
is performed in Section [5l The stability analysis of the vacuum in the brane world 
models is investigated in Section |6l Dynamical dark energy models are considered 
in Section [7] We discuss and conclude our review in Section [5] 

2. Dynamical systems 

2.1. Linear stability of ODE's. In this section we will recall some basics about 
nonlinear systems of ODE's. The reader familiar with these can skip this and the 
next subsection. Our exposition follows closely [2"2"] . 

Let us consider the following system of ordinary differential equations 

(2.1) §-/«, 

where x is a point in some open set U C JR n (or some differentiable manifold, like 
the torus, for instance), and / : U C M n — > M n is a differentiable function. The 
function / is also called a vector field because it assigns a vector f(x) to each point 
in U. If / is a linear function, then (12 . L[) is called a linear system of ODE's. We 
are however concerned with nonlinear systems of ODE's that appear in modeling 
real life phenomena. 

Given the differential equation (|2.ip . let y'(xo) be the solution x(t) with the 
initial condition Xo at t = 0, i.e. 

(2-2) Ax ) = x , |^(x ) = jV(x )), 

at all t for which the solution is defined. It is customary to write ip(t, xo) for y>* (xo). 
The function <^'(xo) is called the flow of the differential equation, and the function 
/ defining the differential equation is called the vector field which generates the 
flow. 

A point x is called a fixed point (or a stationary point) for the flow (p if ip (x) = x. 
The point x it is also called an equilibrium, or a steady-state, or singular point. If 
the flow is obtained as the solution of the differential equation x = /(x), then a 
fixed point is a point for which /(x) = 0. 

A point x is called a periodic point provided there is a time T > such that 
<p T (x) = x and y*(x) ^ x for < t < T. The time T > which satisfies the above 
conditions is called the period of the flow. 

The set O(x) = {<^*(x) : t £ I x } is called the orbit of the point x, where 
I x = (t-,t + ) is the maximal interval of definition of the flow. The orbit of a 
periodic point x is called a periodic orbit. Periodic orbits are also closed orbits 
since the set of points on the orbit is a closed curve. Furthermore, one can define 
the positive semiorbit and the negative semiorbit as + (x) = {<p*(x) : t > 0} and 
0~(x) = {<p*(x) : t < 0}, respectively. 

The positive limit set and the negative limit set of a point x are defined as 

(2.3) A + (x) = {y : there exists a sequence t n — > oo such that lim (/?(£„, x) = y} 
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and 
(2.4) 

A (x) = {y : there exists a sequence t n — > — oo such that lim </j(t„,x) = y}, 

n— >oo 

respectively. 

Obviously, if the forward orbit + (x) is bounded, then A + (x) is connected. 
Let us consider an example from |22j : 

(2.5) x = -y + fj,x(l - x 2 - y 2 ), y = x + (j,y(l - x 2 - y 2 ), 
for \i > 0. In polar coordinates the system becomes 

(2.6) = 1, f = /ir(l-r 2 ). 
The derivative of r satisfies 

(2.7) r > for < r < 1, and f < for 1 < r. 

Therefore, for some initial conditions x ^ (0,0), ro 7^ 0, the trajectory has the 
positive limit set A+(x) = S . Thus solutions forward in time converge onto the 
unique periodic orbit r = 1. 

Such an orbit with trajectories spiraling toward it is called a limit 
cycle. 

For nonlinear systems, one can linearly approximate them near hyperbolic fixed 
points, by means of the Jacobian, such that a linear stability analysis holds good. 

For concreteness, consider a nonlinear differential equation of the type (|2.1|) . and 
let p be a fixed point, namely <p (p) = p for all t. In order to linearize the flow tp* 
near p, let A = (Df)i p = (dfi/dxj)\ p be the Jacobian matrix of / evaluated at p. 
Then the linearized differential equation at p is x = Ax. 

The fixed point p is called hyperbolic if the real part of the eigenvalues of the 
matrix A are nonzero. A hyperbolic fixed point is a sink, or said to be attracting, 
provided the real parts of all eigenvalues of A are negative. A hyperbolic fixed point 
is a source, or said to be repelling, if the real parts of all eigenvalues are positive. 
Finally, a hyperbolic fixed point is called a saddle if it is neither a sink nor a source; 
this would be the case if there exist eigenvalues A+ and A_ such that i?e(A+) > 
and i?e(A_) < 0. A nonhyperbolic fixed point is said to have a center. 

One is often interested in the set of all points whose limit set is a fixed point. 
The stable manifold, or the basin of attraction, of a fixed point p is the set 

(2.8) W s {p) = {x: A+(x)=x}, 
and the unstable manifold of a fixed point p is the set 

(2.9) W u (p) = {x: A-(x)=x}. 

The linearization method of Hartman says that the differential system (|2.ip and 
the linearized system are locally topologically equivalent. Namely, 

Theorem 2.1 ((Hartman) 22J). Consider the differential equation (|2.ip with f 
smooth in an open subset U of M n containing the origin. If the origin x = is 
a hyperbolic fixed point, then in a small neighborhood of x = 0, there exist stable 
and unstable manifolds W s and W u with the same dimensions as the stable and 
unstable manifolds E s and E u of the linearized system, respectively, where W s and 
W u are tangent to E s and E u at x = 0. 

If the fixed point is not hyperbolic, then other methods have to be used. 
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2.2. Planar Nonlinear differential equations systems. The phase-portrait of 
a two-dimensional differential equation x = /(x) is a drawing of the solution curves 
with the direction of increasing time indicated. In an abstract sense, the phase 
portrait is the drawing of all solutions, but in practice it only includes the repre- 
sentative trajectories. The phase space is the domain of all x's considered. 

Definition. Let p be a fixed point of the two-dimensional system x = /(x), and 
denote by Ai, A2 the two eigenvalues of A := (Df)\ p . The following classification 
of the fixed point p is standard. 

(1) Ai, A2 are real and distinct. 

1.1. Ai • A2 > (the eigenvalues have the same sign): p is called a 
node or type I singularity; that is, every orbit tends to the origin in a 
definite direction as t — > 00. 

1.1.1. Ai, A2 > : p is an unstable node. 

1.1.2. Ai, A2 < : p is a stable node. 

1.2. Ai • A2 < (the eigenvalues have different signs): p is an unstable 
fixed point, or a saddle point singularity 

(2) Ai, A2 are complex, i.e. Ai,2 = ct ± ij3, (3^0. 

2.1. a / 0: p is a spiral, or a focus, that is, the solutions approach the 
origin as t — > 00, but not from a definite direction. 

2.1.1. a < 0: p is a stable focus. 

2.1.2. a > 0: p is an unstable focus. 

2.2. a — 0: p is a center, that means it is not stable in the usual sense, and 
we have to look at higher order derivatives. 

(3) Ai, A2 are equal, i.e. Ai = A2 = A. 

3.1. If there are two linearly independent eigenvectors, we have a star singu- 
larity, or a stable singular node (these are simple straight lines through 
the origin). 

3.2. If there is only one linearly independent eigenvector, we have an im- 
proper node, or unstable degenerate node. 
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FIGURE 1. Stable (left) versus unstable (right) nodes. 
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Figure 2. Stable (left) versus unstable (right) spirals. 





Figure 3. Saddle point singularity (left) and a center (right). 



Remark 1. As one can easily remark from the above classification, linear stability 
means that the fixed point p is an attractor for the nearby solutions, or, geomet- 
rically, that the trajectories are "converging". For the sake of simplicity, let us 
consider here a one dimensional autonomous differential equation x — f(x) in the 
Euclidean plane (R, || • ||) with a fixed point p. This means that if a solution starts 
at this point, it remains there forever, i.e. p(t) = p, for all t. 

Strictly speaking, the critical point p is called stable if given e > 0, there is a 
5 > such that for all t > t , \ \x(t)—p(t)\ \ < e whenever t > t , \ \x(t )—p(to)\ \ < S, 
where x(t) is a solution of x = f(x). One can see that this definition is equivalent 
to the fact that f'{x ) < 0, as predicted by the general case. Indeed, let us consider 
a small perturbation £ (t) away from the critical point p, i.e. 



(2.10) 
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Figure 4. Stable singular node (left) versus unstable degenerate 
node (right). 



We will analyze now the tendency of the perturbation £(t) to grow or decay in time. 
By taking the derivative of (|2.10[) it follows 

(2.11) £(*)=i(*) = /(*) = /(p + 0> 

and by Taylor expansion, we get 

(2-12) J(t)=/(P)+J/'(P) + Y/"(P) + -- 

We are interested in linear analysis, therefore the higher order terms can be 
neglected. Taking into account that f(p) = 0, we obtain 

(2.13) i = tf'{p) 

We can immediately conclude from here that 

(1) if f'ip) < 0, then the perturbation £(t) decays exponentially, i.e. trajecto- 
ries are converging, and 

(2) if f'(p) > 0, then the perturbation £(t) grows exponentially, i.e. trajectories 
are diverging. 

The higher dimensional case can be analogously treated leading to the classifi- 
cation above. 

2.3. Limit Cycles. Roughly speaking, a limit cycle is an isolated periodic solution. 

It is known that for a planar ODE system, the trajectories tend to a critical 
point, a closed orbit, or to infinity. 

In order to discuss the stability of limit cycles, let us recall the notion of invariant 
set. A subset S £ R 2 is called invariant with respect to the flow if, for any x £ S, 
it follows <^*(x) £ S, for all t £ I x . Similarly, a subset S £ M 2 is called positively 
invariant and negatively invariant with respect to the flow ip t if, for any x £ S, it 
follows + (x) £ S, and 0~(x) £ S, respectively. 

Then, a limit cycle, say $, is called 
(1) a stable limit cycle if A + (x) = $, for all x in some neighborhood, this 
obviously means that the nearby trajectories are attracted to the limit cycle; 
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(2) an unstable limit cycle if A - (x) = <!>, for all x in some neighborhood, 
meaning that the nearby trajectories are repelled from the limit cycle; 

(3) a semistable limit cycle if it is attracting on one side and repelling on the 
other. 

Recall here an useful tool for the study of existence of limit cycles in plane, 
namely the Poincare-Bendixon Theorem. 

Theorem 2.2 (Poincare-Bendixon). Suppose (9 + x is contained in a bounded region 
in which there are only finitely many critical points. Then A + (0 + ) is one of the 
following 

(1) a single critical point; 

(2) a single closed orbit; 

(3) a graphic, i.e. critical points joined by heteroclinic orbits. 

In concrete problems, the following corollary is usually used 

Corollary 2.3. Let S be a bounded closed set containing no critical points and 
suppose that S is positively invariant. Then, there exists a limit cycle entirely 
contained in S. 

The eigenvalues of an appropriate matrix can be used to determine the stability 
of periodic orbits, or limit cycles. This is the so-called Floquet theory. 

We are interested in the study of dynamics of a flow in a neighborhood of a 
periodic orbit. One method to do this is to follow the nearby trajectories as they 
make one circuit around the periodic orbit. 

Let tp be a periodic orbit of period T and p € tp, i.e. f T {p) = p. Then for some 
k € N, the fc-th coordinate function of the vector field must not vanish at p, i.e. 
there exists a k such that fk(p) ^ 0. We consider the hyperplanc S = {x| Xfe = pk}, 
called a cross section or transversal at p. 

For x £ £ near p the flow tp *(x) returns to £ in time t(x), which is roughly 
equal to T (this follows from the Implicit Function Theorem). 

Let V C £ be an open set on which t(x) is a diflerentiable function. The first 
return map or Poincare map V : £ ->• £ is defined to be V(k) = tp T W(x) for x G V, 
see Fig. [5] 

It can be seen that the Poincare map is diflerentiable. If ip is a periodic orbit of 
period T with p £ tp then it can be shown that the eigenvalues of (D ( p T )\ x =p are 
1, Ax, A n _i (the eigenvector of 1 is f(p))- The n — 1 eigenvalues Ai, A n _i are 
called characteristic multipliers of the periodic orbit ip. It is remarkable that these 
depend only on tp but not on the point p. 

The central result of Floquet theory says that the eigenvalues of the Poincare 
map's Jacobian matrix at p (i.e. (DV)\ x =p) coincide with the characteristic multi- 
pliers of the limit cycle. Based on this result, the stability of limit cycles has been 
formulated in terms of their characteristic multipliers. 

A limit cycle tp is hyperbolic provided \Xj\ ^ 1 for all characteristic multipliers 
(1 < j < n — 1). It is a hyperbolic stable limit cycle, a periodic attractor, or a 
periodic sink, if all \Xj\ < 1. It is a hyperbolic unstable stable limit cycle, a periodic 
repeller, or a periodic source, if all |A,| > 1. A hyperbolic periodic orbit which is 
neither a source or a sink is called a saddle limit cycle. 

Remark 2 (The planar systems case). For a two dimensional differential equation 
the Poincare map can sometimes be explicitly calculated. Indeed, let us consider 
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Figure 5. The Poincare map. 



the differential equation 

(2.14) x = X(x,y) y = Y(x,y) 

and 



(2.15) V(x,y) 



X(x,y) 
Y(x,y) 



Let us consider the transversal E and let E' C E be the open set where Poincare 
map is defined, V : E' — > E. Since we are in a plane, E is a curve which can be 
reparametrized by 7 : / — > E with j(I') = E', I' C / and |7(s)| = 1. As in the 
general case, let r(q) be the return time for q G E', so V(q) = ip T ^(q). 

If we assume V(qo) — qo and j(so) — qo, then 

(2.16) (V o 7 )'( fl0 ) = exp[ / (<?0) (^ y)o^(g )di] 

Jo 

(see [12] for details). 

In this case, one can define the characteristic multiplier M to be 

dV 

2.17 M := — , 

where g* is a fixed point of the Poincare map V corresponding to a limit cycle, say 
If \M\ < 1, then $ is a hyperbolic stable limit cycle, if \M\ > 1, then $ is an 
unstable limit cycle, and if \M\ = 1, 4-? 7^ 0, then $ is a semistable limit cycle. 



Let us also recall here that, for an ODE system depending on a parameter k, 



i.e. 



d-x 

(2-18) M =f{x ' k) > 

a value, say is called bifurcation value if at the point (x, k) the stability of the 
solution changes. For planar systems, there are several types of bifurcation at 
nonhyperbolic critical points. A bifurcation concerning the stability of a limit cycle 
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is called a Hopf bifurcation. This is the kind of bifurcation we are most interested 
in. Basically, there are two types of Hopf bifurcation: 

(1) the supercritical Hopf bifurcation when stable limit cycles are created about 
an unstable critical point, and 

(2) the subcritical Hopf bifurcation when an unstable limit cycle is created about 
a stable critical point. 

Finally, we mention that the problem of determining the maximum number of 
limit cycles for planar differential systems is still an open problem. This is the 
second part of the 23 rd Hilbert problem and only little progress has been made to 
date. 



3. Kosambi-Cartan-Chern (KCC) theory and Jacobi stability 

3.1. General Theory. We recall the basics of KCC-theory to be used in the sequel. 
Our exposition follows 23 . 

Let M be a real, smooth n-dimensional manifold and let TM be its tangent 
bundle (in the most cases in the present exposition we will consider the case M — 
M. n ). Let u — (x,y) be a point in TM, where x = (x 1 , x 2 , x n ), and y = 

(y\y 2 f). 

A time independent coordinate change on TM is given by 

(3.1) Q^.i i E {1, 2, n} . 

Let us also recall that the equations of motion of a Lagrangian mechanical system 
(M,L,Fi) can be described by the famous Euler- Lagrange equations 

d dL dL _ 

(3.2) dtW-^- Fl ^l, 2 ,...,n, 

where L is the Lagrangian of M. , and Fi are the external forces [19l [18] . 

For a regular Lagrangian L, the Euler-Lagrange equations given by Eq. (13. 2[) 
are equivalent to a system of second-order differential equations 

(3.3) ^! + 2G 4 (x, 2/ ) = 0, ie{l,2,..,n}, 

where G z (x, y) are smooth functions defined in a local system of coordinates on 
TM. 

The system given by Eq. (|3.3j) is equivalent to a vector field (semispray) S on 
TM, where 

(3.4) S = y l —~2G l ix\y\t)— -, 

which determines a non-linear connection N on TM with local coefficients TVj 
defined by [18] 
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More generally, one can start from an arbitrary system of second-order differen- 
tial equations on the form (|3.3p . where no a priori given Lagrangean function is 
assumed, and study the behavior of its trajectories by analogy with the trajectories 
of the Euler-Lagrange equations. 

For a non-singular coordinate transformations given by Eq. (|3.1I) . we define 
the KCC-covariant differential of a vector field £ = £*(;c)g§r in an open subset 

n C R n x R n as [TJ d [23 HO 
(3.6) E1L = ^L + N ia. 



dt dt 



For £ 1 = y l we obtain 



(3.7) ^f- = N^y 1 - 2G l = -e\ 

dt J 

The contravariant vector field e 1 on VL is called the first KCC invariant. We point 
out that the term 'invariant' refers to the fact that e l has geometrical character, 
namely, with respect to a coordinate transformation (13.11) . we have 

< 3 - 8 > ? = jS* 

Let us now vary the trajectories x l (t) of the system Q3.3P into nearby ones ac- 
cording to 

(3.9) x i (t)=x i (t) + r ] e(t), 

where \rj\ is a small parameter and are the components of some contravariant 
vector field defined along the path x z (t). Substituting Eqs. Q3.9p into Eqs. Q3.3P 
and taking the limit 77 — > we obtain the variational equations [J [H [531 HI] 

By using the KCC-covariant differential we can write Eq. (13.101) in the covariant 
form 

D 2 P ■ ■ 

(3.11) *r = W 



where we have denoted 



dG i , ■ ,dN) . , 



(3.12) P] = -2— - 2G'G}; + y'^f + JV?JV^ 

and Gjj = dNj/dy 1 is called the Berwald connection [TJ [H US [M] . Eq. (|3~TT|) 
is called the Jacobi equations, or i/ie variation equations attached to the system 
of second order differential equations, and P? is called i/ie second KCC-invariant 
or i/ie deviation curvature tensor. When the system (|3.3p describes the geodesic 
equations in either Riemann or Finsler geometry, Eq. p. lip is the usual Jacobi 
equation. 

The third, fourth and fifth invariants of the system (|3.3p are given by [2] 



(313) p* =N^i Ml p< = d Ik D i = d 3i 

The third invariant is interpreted as a torsion tensor, while the fourth and fifth 
invariants are the Riemann- Christoffel curvature tensor, and the Douglas tensor, 
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respectively 2 . These tensors always exist and they describe the geometrical prop- 
erties of a system of second-order differential equations (|3.3|) ( [2] 123] [24] ) . 

Definition. The trajectories of (|3.3[) are Jacobi stable if and only if the real parts 
of the eigenvalues of the deviation tensor Pj are strictly negative everywhere, and 
Jacobi unstable, otherwise. 

Let us discuss the geometrical meaning of this definition in the context of an 
Euclidean, Riemannian or Finslerian structure. 

Since in many physical applications we are interested in the behavior of the 
trajectories of the system (|3.3j) in a vicinity of a point x l (to), where for simplicity 
one can take to = 0, we will consider the trajectories x l = x l (t) as curves in the 
Euclidean space (R n , (., .)), where (., .) is the canonical inner product of R n . As for 
the deviation vector £ we assume that it satisfies the initial conditions f (0) = O 
and i (0) = W ^ O, where O € R n is the null vector [251121] . 

For any two vectors X,Y e R n we define an adapted inner product ((., .)) to 
the deviation tensor £ by ({X,Y}) := 1/ (W, W) ■ (X,Y). We also have ||VF|| 2 := 
{{W,W)) = 1. 

Thus, the focusing tendency of the trajectories around to = can be described 
as follows: if ||£ (t)\\ < t 2 , t w + , the trajectories are bunching together, while 
if > t 2 , t ^ + , the trajectories are dispersing (23j [24]. In terms of the 

deviation curvature tensor the focusing tendency of the trajectories can be described 
as follows: The trajectories of the system of equations p.3[) arc bunching together 
for t s» + if and only if the real part of the eigenvalues of (0) are strictly 
negative, and they are dispersing if and only if the real part of the eigenvalues of 
P] (0) are strict positive p3] [24]. 

The focussing behavior of the trajectories near the origin is represented in Fig. [6] 



~x\t) ¥{t) 




Figure 6. Behavior of trajectories near zero. 

In the more general case of a typical SODE of the form 13.31 (not necessarily 
associated to any metric), let us remark that Jacobi instability via equation (|3.Iip 
means exponential growth of the deviation vector field £, while Jacobi stability 
means that £ have an oscillatory variation expressed by a linear combination of 
trigonometric functions sin and cos. 
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3.2. The relation between linear stability and Jacobi stability. It would 
be interesting to correlate linear stability with Jacobi stability. In other words, to 
compare the signs of the eigenvalues of the Jacobian matrix J at a fixed point with 
the signs of the eigenvalues of the deviation curvature tensor Pf evaluated at the 
same point. Even though this should be possible in the general case, we give here 
only a discussion for the 2-dimensional case. 
Let us consider the following system of ODE: 

. , du ,. N dv . . 

(3-14) — = f( u ,v) —=g(u,v) 

such that the point (0, 0) is a fixed point, i.e. /(0, 0) = g(0, 0) = 0. In general, even 
though the fixed point is (ito, vq), by a change of variables u — u — uq, v = v — vq, 
one can obtain in most of the cases (0,0) as a fixed point. We denote by J the 
Jacobian matrix of (I3.14p . i.e. 



(3.15) J(u,v) 



fu fv 

9u 9v 



where the subscripts indicate partial derivatives with respect to u and v. 
The characteristic equation reads 

(3.16) A 2 - trA + det A = 0, 

where tr A and det A are the trace and the determinant of the matrix A := J 

(0,0)' 

respectively. 

The signs of the trace, determinant of A, and of the discriminant A = (f u — 
9v) 2 + ^fv9u — {trA) 2 — 4 det A give the linear stability of the fixed point (0, 0) as 
described in the previous section. 

Let us show how is possible to apply the KCC theory to the ODE system (|3.14l) . 
By elimination of one of the variables, we can transform (13.141) into a SODE of 
the form (13.31) and compute the deviation curvature. It does not matter which 
of the two variables we eliminate, because the resulting differential equations are 
topologically conjugate ([22]). 

Let us relabel v as x, and g{u,v) as y, and let us assume that g u \(a.o) ^ 0. The 
variable to be eliminated, u in this case, should be chosen such that this condition 
holds. Since (u, v) — (0, 0) is a fixed point, the Theorem of Implicit functions 
implies that the equation g(u, x) — y — has a solution u — u(x, y) in the vicinity 
of (x, y) — (0,0). Now, 'x = g = g u f + g v y. Hence we obtain an autonomous 
one-dimensional case of the SODE (|3.3I) . namely x 1 + g 1 = 0, where 

(3.17) g x (x, y) = -g u (u(x, y), x) f(u(x, y), x) - g v (u(x, y), x) y. 

Evaluating now the curvature tensor Pj 1 from (j3. 1 1[) at the fixed point (0,0), we 
obtain after some computation: 

Theorem 3.1. Let us consider the ODE (|3.14l) with the fixed point p—(0,0) such 
that g«|( ,o) 7^ 0. 
Then, 

(3.18) 4P 1 1 | (0 ,o) = -4^11(0,0) + (ff;i) 2 |(o,o) = A := {trA) 2 - 4 det A. 
In this case, the trajectory v = v(t) is Jacobi stable if and only if A < 0. 
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It can be checked that the second equality would not change if one had used 
the first equation of (|3.14|) to eliminate the variable v instead by labeling u as x, 
provided /„|( , ) 7^ 0. 

In other words, we can conclude that in the case g u \(o.o) 0, fv\(o,o) both 
trajectories u = u(t) and v — v(t) are Jacobi stable if and only if A < 0. 

Now, straightforward calculations give the following 

Corollary 3.2. Let us consider the ODE (|3.14j) with the fixed point p=(0,0) such 
that g u \(o,o) 0. 

Then, the Jacobian J estimated at p has complex eigenvalues if and only if p is 
a Jacobi stable point. 




Figure 7. Linear stability. 

One can represent the linear stability as in the Fig. and by combining with 
the definition of Jacobi stability, in the conditions of Theorem above, we obtain the 
following relations 



(1) Region A. 



(2) Region B. 



(3) Region C. 



A > Jacobi unstable 

S = trA > Unstable node 

P = dctA > 



A < Jacobi stable 

S = trA > Unstable focus 

P = dctA > 



A < Jacobi stable 

S = trA < Stable focus 

P = detA > 
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(4) Region D. 



(5) Region E. 



A > Jacobi unstable 

S = tiA < Stable node 

P = Act A > 



A > Jacobi unstable 

P = dctA < Saddle point 

Let us point out that the stability in Jacobi sense refers to a linear stability type 
of the trajectories in the curved space endowed with a nonlinear connection and 
a curvature tensor, as described above. Here the role of usual partial derivative is 
played by the covariant derivative along the flow. This obviously leads to difference 
in the meaning of linear stability and Jacobi stability. 

3.3. Example: The Brusselator. We present a first example of a dynamical 
system that has a limit cycle behavior, but which is not Jacobi stable in its entirety. 
In other words, the Brusselator is not robust in its entirety. 

The Brusselator is a simple biological system proposed by Prigogine and Lefever 
in 1968 [21] (the name Brusselator came from Brussels, the city where they lived). 
The chemical reactions are: 



(3.19) 



where the fc's are rate constants, and the reactant concentrations of A and B 
are kept constant. Then, the law of mass action leads to the following system of 
differential equations: 

du 





A 


fci 


x, 


BA 


-X 


k 2 


Y + D 


2X - 


fr 


fc 3 


3X, 




X 


k>4 


E 



(3.20) 



1 — (b + l)u + au 2 v —: f(u, v) 
at 



—— = bu — au 2 v =: q(u, v), 
dt Jy ' 



where u and v correspond to the concentrations of X and Y, respectively, and a, 
and b arc positive constants. 

3.3.1. Stability of fixed points. From the equations (|3. 191) it follows that there is a 

unique fixed point S(uq = l,Vo = —)• 

a 

Simple computations show that the Jacobian evaluated at the fixed point S is 
given by 

= J\(u ,v ) 

— (b + 1) + 2auv au \, 
b - 2auv -au 2 ) l«o=l,^o=-| 

6-1 a 

—b —a 
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Region 


A 


B 


C 


D 


Linear 
stability 


Unstable 
node 


Unstable 
focus 


Stable 
focus 


Stable 
node 


Jacobi 
stability 


Jacobi 
unstable 


Jacobi 
stable 


Jacobi 
stable 


Jacobi 
unstable 



Table 1 . Linear stability versus Jacobi stability for the Brussela- 
tor system. 



and hence we have 

(3.21) ttA = b-l-a, detA = a, A = (b - 1 - a) 2 - 4a. 




1 2 3 4 5 6 7 



Figure 8. The parametric regions indicating the linear stability 
domains of the fixed point 5 of the Brusselator. 

A graphical representation of these functions of a and b is given in Fig. [3] Taking 
into account that detA = a > 0, a similar analysis with the one done in the previous 
section gives the following table. 

We would like to illustrate this discussion with a concrete motivation of the 
apparent discrepancies between linear and Jacobi stability of fixed points. For in- 
stance, if the fixed points 5 belongs to Region D, then it must be a stable node, 
i.e. trajectories are attracted to 5 in the phase plane, while they are divergent 
in the Jacobi stability space. A quick look at the previous section show the main 
differences between the variation of the deviation vector field (t) given by the for- 
mulas (|3.10[> and (j3.11|) . respectively. Obviously the local behavior of the solutions 
of these two equations are completely different. See figures 

The same kind of analysis can show how the deviation vector field variate in all 
four different regions explaining in this way the Table I given above. 

3.3.2. Limit cycles. The following result is well known about the Brusselator (see 
for example [20] ) 



Corollary 3.3. The fixed point 5(1, b/a) is 
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Figure 9. Variation of the deviation vector field given by 
(pUP) (left) and (jgXEJ) (right), for a = 4, b = 0.5. 



o an unstable node if and only if b > (*/a + l) 2 , and 

o an unstable spiral if and only if b < (\fa + l) 2 . 

Moreover, the point b c = a+1 is a critical point in the sense of Hop f, and for b > b c 
the system exhibits a limit cycle behavior. 

Indeed, it can be shown that in this case there is a confined closed set that con- 
tains the fixed point 5(1, — ) in its interior. From the Poincare-Bendixon theorem, 
it follows that the system exhibits a limit cycle. 

We have seen that if the fixed point 5(1, b/a) belongs to one of the regions A or 
B, then it must behave as a limit cycle. 

From Table I one can easily see that the Jacobi stability of the limit cycle changes 
when moving from the region A to B. 

The usual bifurcation analysis does not give any information about the stability 
of the periodic solutions. We are going to estimate the stability of such a solution 
by means of Jacobi stability. 

First, we remark that because of trA > 0, and detA > 0, from Corollary 3.5 it 
follows: 

Proposition 4.2 

• If the fixed point 5(1, ^) is an unstable spiral, then it is in the Jacobi stabil- 
ity region. Conversely, if it is in the Jacobi stability region, and b > a+1, 
then it is an unstable spiral. 

• If the fixed point 5(1, — ) is an unstable node, then it is in the Jacobi in- 
stability region. Conversely, if it is in the Jacobi instability region, and 
b > a + I, then it is an unstable node. 

This can also be seen directly from 4P 1 1 | ( - 1 ^ = A = (trA) 2 — AdetA. 

Therefore, in an open region around 5, the transient trajectories are Jacobi 
stable, so we can conclude that in this region the system is robust and it becomes 
fragile outside it. This conclusion is in accord with the usual stability analysis done 
with Floquet theory. 
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The next sections introduce various physically motivated dynamical systems to 
which the previously discussed techniques can be applied. While following the linear 
stability results are mainly know, most of the remaining analysis is new. 

4. Newtonian polytropes - Lane-Emden equation 

It is well known that the Lane-Emden or Emden-Fowler equation can be rewrit- 
ten as an autonomous system of equations by introducing a new set of variables. 
The Lane-Emden equations appear very naturally in the study os compact astro- 
physical objects in Newtonian gravity. 

4.1. The Lane-Emden equation and linear stability. The properties of the 
stable Newtonian compact astrophysical objects can be completely described by 
the gravitational structure equations, which are given by: 

(4.1) — ±-L = 4:T T pr A , 

dr 

dp Gm(r) 

( 4 -2) -j- = — P, 

ar r z 

where p > is the density, p > is the thermodynamical pressure, m(r) denotes 
the mass inside radius r, satisfying the condition m(r) > 0, V?' > 0. By eliminating 
the mass between Eqs. (|4.1|) and (|4.2|) we obtain the equation [14] 

1 d fr 2 dp^ 
dr \ p dr / 

To close this structure equation, one assumed an equation of state to specify the 
matter content, p — p{p). We assume that the equation of state is continuous and it 
is sufficiently smooth for all p > 0. A physically meaningful solution of Eqs. (|4.1[) - 
(I4.3|) is possible only when boundary conditions are imposed. We require that 
the interior of any matter distribution be free of singularities, which imposes the 
condition m(r) — > as r — > 0. At the center of the star the other boundary 
conditions for Eqs. (|4.1[) - (|4.3[) are p(0) = p c , p(0) — p c , where p c and p c are the 
central density and pressure, respectively. The radius R of the star is determined 
by the vanishing pressure surface p(R) = 0. 

An isotropic bounded fluid distribution for which the pressure and the density 
are related by a power law of the form 

(4.4) p = Kp 1+1 / n , 

where K and n are constants, plays an important role in astrophysics, n > is 
called the polytropic index. It is convenient to represent the density in terms of a 
new dimensionless variable 9 defined by 

(4.5) P = p c n , 

giving for the pressure p = K p\ +1 ^ n 9 n+1 . By rescaling the radial coordinate as 

I 1/2 



(4-3) ir,iz[--t.)=-^G P . 



(n + 1) Kp X J n 1 /4nG , ±oo, for polytropic stars Eq. (I4.3[) 



takes the form of the Lane-Emden equation of index n, 
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The initial conditions for the Lane-Emden equation are 9(0) = 1 and 0'(0) = 0, 
respectively, where the prime denotes the derivative with respect to the variable £. 
By introducing a set of new variables (w, t) defined as [13] 

(4.7) 9(0 = Be /{1 - n) w(0, £ = e-*, 

where B > is a constant, the Lane-Emden equation is transformed into the 
following second order differential equation, 

drw 5 — n dw 2 (3 — n) 
~dlF + n-l~dt 



(4.8) 



-w + S n_1 iy n = 0, n^l,±oo. 



(n-iy 

In the following we will restrict the range of the polytropic index n to the range 
n G (l,+oo). From a mathematical point of view Eq. (14.8[) is a second order 
non-linear differential equation of the form d 2 w/dt 2 + 2G 1 (w,dw/dt) = 0, with 



(4.9) 



G 1 f w, 



dw\ 
It J 



5 — n dw 2 (3 — n) 
n — 1 dt 



w + B n - 1 w n 



(4.10) 



= -2G\w,q). 



(n-iy 

The function G 1 has the properties G x (0,0) = and G 1 (w,0) ^ for w ^ 0, 
respectively. 

In order to study the stability of the equilibrium points of Eq. (|4.8I) , we rewrite 
it in the form of an autonomous dynamical system, by introducing a new variable 
q = dw/dt, so that the Lane-Emden equation is equivalent to the following first 
order autonomous system of equations 

dw dq 
~dt^ q ' dt 

The critical points of this dynamical system are given by q = 0, and by the solutions 
of the equation 

(4.11) G 1 (w,0)=Q. 

Therefore for the critical points X n = (wo,qo) of the system ()4.10|) we find the 
values 

X o = (0,0), ne(l,+oo), 

-il/Cn-l) 

•Y„ - i - — 4 ,0 | , t»€(3,oo), 



(4.12) 
(4.13) 



(n-iy 



(4.14) 



, : 2/(n-l) 



2 (n - 3) 
(n-lf 



l/(n-l) 



.0 



n e (1,3). 



4.2. Jacobi stability analysis of Newtonian polytropic fluid spheres. Fol- 
lowing our previous work |J we denote x := w and y := dx/dt = dw/dt. Then 
Eqs. (|4.8[) can be written as 



(4.15) 
where 

(4.16) 



G 1 (x,y) 



^+2G' (x,y) = 0, 



n 2 (3 - n) 
-rV+ — Jrx + B h 



(n-iy 



(x) n 
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respectively, which can now be studied by means of the KCC theory. The follow- 
ing treatment can also be found in one of our previous works [4]. As a first step 
in the KCC stability analysis of the Newtonian polytropes we obtain the nonlin- 
ear connections = dG 1 /dy, associated to Eqs. f|4.15[) . and which is given by 
Nl(x, y) = (5 — n) /2(n — 1). For the Lane-Emden equation the associated Berwald 
connection G\ 1 = dNl/dy is identically equal to zero G\ x = 0. Finally, the second 
KCC invariant, or the deviation curvature tensor Pi, defined as 



(4.17) 

is given by 

(4.18) 



9 N 1 »,1 »rl dN l 



Pl(x,v) = - A 



dx 

nB n - 1 (x) n - 1 



dt 



or, in the initial variables, by 
(4.19) Pi(w,q) 



nB 



™- 1 w; rl - 1 . 



Evaluating Pl(w,q) at the critical points X n , given by Eqs. (|4.12[) (|4.13[) . we 
obtain 

P 1 1 (X o ) = i>0, ne(l,oo), 



(4.20) 
(4.21) 

(4.22) 



Pi (Xi n ) — 



Pi (Xn) = 



9n 2 - 26n + 1 
4(n-l) 2 : 

-7n 2 + Tin + 1 



n E (1,3). 



n 6 (3, oo). 



4(n-l) 2 ' 

The behavior of the functions Pj (n) is represented in Fig [TUJ 




1.25 1.5 1.75 2 2.25 2.5 2.75 




(4.23) 



Figure 10. Left: The deviation curvature tensor Pj 1 (Xj n ) as a 
function of n € (1 5 3). Right: The deviation curvature tensor 
Pi (X n ) as a function of n > 3. 

In the initial variables the deviation curvature tensor Pi is given by 

1 



Pl&e) = ~ n ee n -\ 



Pi can be expressed in a simple form with the use of Milne's homological variables 
(u,v), defined as [151 [Tl] 
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with the use of whom we obtain 

(4.25) Pi(u,v) = - -nuv. 

From a physical point of view u(r), defined as 

(4.26) u(r) = dlnm(r)/dlnr = 3p(r)/p(r), 

is equal to three times the density of the star at point r divided by the mean density 
of matter p(r) contained in the sphere of radius r. The variable v{r), defined as 
v(r) = -d1np(r)/dlnr = (3/2) [Gm(r)/r] / [3p(r)/2p(r)], is 3/2 of the ratio of the 
absolute value of the potential energy \E g \ = Gm(r)/r and the internal energy per 
unit mass Ei = (3/2)p/p, so that v(r) = (3/2) \E g \ /Ei [H]. In terms of these 
physical variables the deviation curvature tensor is given by 

(4.27) Pl{r)- 1 3n "WI^I 



4 2 p(r) Ei 

The condition of the Jacobi stability of the trajectories of the dynamical system 
describing a Newtonian polytropic star, Pj 1 < 0, can be formulated as 

(4-28) §- < 6n^. 

\ E g\ P{r) 



5. The general relativistic static fluid sphere 

The properties of an isotropic compact general relativistic object can be com- 
pletely described by the gravitational structure equations, which are given by [8]: 

(5.1) — =4npr , 

dr 



dp _ (p + p) (to + 4irr 3 p) 
( ' dr ~ r 2 (1 - 2f ) 

where m(r) is the mass inside radius r satisfying the condition m(r) > 0,Vr > 0. 
To close the field equations an equation of state p = p(p) must also be given. We 
assume that the barotropic equation of state is continuous and it is sufficiently 
smooth for p > 0. In the Newtonian limit Eq. (|5.2[) reduces to the expression given 
by Eq. 

A physically meaningful solution of Eqs. (|5.1j) ~ (|5.2j) is only possible when bound- 
ary conditions are imposed. We require that the interior of any matter distribution 
be free of singularities, which imposes the condition m(r) — > as r — > 0. At the 
center of the star the other boundary conditions for Eqs. (|5.1I) - (|5.2I) are 

(5.3) p(P)=Pc, p(0)=Pc 

where p c and p c are the central density and pressure, respectively. The radius R of 
the star is determined by the boundary condition p(R) — 0. Depending on the form 
of the equation of state of the matter several fluid star models can be considered. 
By introducing the set of the homologous variables P, p, M and 9, defined as 



(5.4) tt(v) = Airr 2 p(p(v)), v — Aitr 2 p, u — m/r, t = bx(r/R), 
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the structure equations of the star take the form 

(5.5) J = 

^ 2tt - ( t; + 7r ( t; ))( M + 7r ( 7; )) 

1 ' J dt ^ 1 - 2m 

On the boundary of the star, where r = R we have t — 0. At the center of the star, 
where r — > 0, we have f — >• oo. In order to close the system of equations (|5.5|) ~ (|5.6|) . 
the equation of state of the matter must be given in the form n = tt(v). 

5.1. Compact objects and linear stability. As we have seen in Section \'3.1l 
the stability properties of a dynamical system of second order differential equations 
can be investigated by using the KCC geometric approach. The KCC theory gives 
the possibility of systematically investigating the Jacobi stability properties of the 
general relativistic compact objects in a rigorous framework. As a first example of 
the application of the KCC theory to the fluid sphere model we consider the case 
of in which the equation of state of the fluid is given by a linear relation between 
the pressure and the energy density, so that 

(5.7) p=( 1 -l)p, K 7 <2. 

The equation of state of the homologous variables is ir(v) = (7 — l)v, and the 
sound speed in the medium is given by c 2 s =7—1. The structure equations of 
the isotropic star with a linear barotropic equation of state take the form of an 
autonomous system jSJ, given by 

(5.8) f _ „-«, 

/- „s dv „ v / 57 — 4 

(5 - 9) Tt = 2v -—){ 2 -^~^T u 

By eliminating \x between these equations we obtain the second order differential 
equation for M in the form 



(5.10) 



cPu du (u + du/dt) 
~d¥ + ~dt 1 - 2u 



7" + 47 — 4 du 

u — 7— - 

7-1 ' dt 



= 0. 



We give here a short review of the linear analysis of the autonomous system 
given by Eqs. (15.81) and (|5.9[) . respectively. One can find a similar analysis in [8]. 
There are two steady states for this dynamical system, namely 

(5.11) («o,«o) = (0,0), 

2(7-1) 2( 7 -l) 



(5-12) (ui,*!) = v , . 

V 7 Z + 47 — 4 7^+47 

respectively. 

Let us remark two things about these steady states. Firstly, one can see that 
(uo,Vo) is not desirable because of the conditions u(t) > 0, v(t) > imposed on 
the system by l|5.4|) . and therefore the only acceptable steady state is (ui, One 
should also remark that for u = 1/2, the system has a singularity, so actually we are 
concerned only with < u < 1/2. Secondly, for 7 = —2 ± 2V2, the denominator of 
(ui,vi) vanishes. However, taking into account the condition 1 < 7 < 2, we remark 
that (ui,V\) is well defined in our case. 
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The Jacobian matrix of Eqs. (|5 .8[) and (|5.9p is given by 

/ 2 7 (7-1)m-(57-4)M-2( 7 -1) 7 M [1+2( 7 -1)m] \ 

(5.13) J = I (7-l)(2M-l) ( 7 -l)(2M-l)* j ; 

and by evaluating it at (tti,Wi) we obtain 

/ 2(7-1) 2(57-4) \ 

(5.14) A 1 = J(u 1 ,v 1 )= I y J. 

The eigenvalues of the matrix A\ are the solutions of the characteristic equation 

(5.15) A 2 - A-tr(A x )+detAi =0, 
where 

(5.16) ^0=-^, ^ Al = 2 ^ +A r A \ 

7 7 

are the trace and determinant of the matrix A±. Moreover, the discriminant of the 
characteristic equation (|5.15[) is given by 

(5.17) A = (til Ax)) - 4det = ^— — . 

7 

Remark that for 1 < 7 < 2 we always have A < 0, i.e. the eigenvalues of the 
matrix A\ are complex. Moreover, since tr(Ai) < for the given range of 7, it 
follows that the steady state (u\,v\) is a stable spiral, as shown in Fig. [TT] 




Remark that for 1 < 7 < 2 there are no bifurcations points, i.e. the stable spirals 
are very robust. 
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5.2. Jacobi stability analysis. The study of the Jacobi stability analysis of the 
general relativistic fluid sphere with a linear barotropic equation of state is based 
on the direct study of the second order differential equation given by Eq. (|5.10l) . 
In order to simplify the notation we introduce the new variables x := u and y :— 
dx/dt — du/dt. Therefore Eq. (|5.10[) can be written as 



(5.18) 



cPx 
d6 2 



I) 



jx + y) 

l-2x 



2- 



T + 4 7 - 

r 



7- 



= 0. 



In other words, we can express the system of the structure equations of the compact 
general relativistic star with a linear barotropic equation of state, given by Eqs. (|5.8p 
and (|5.9[) . respectively, by the SODE 



(5.19) 
where 
(5.20) 



^±+2G 1 (x,y)=0, 



2G 1 (x,y) = y- 



l-2x 



2- 



T + 47 - 
7-1 



-x — jy 



As a first step in the KCC stability analysis of the general relativistic fluid 
sphere with linear barotropic equation of state we obtain the nonlinear connection 
Nl associated to Eq. (|5.19[) . and which is given by 



(5.21) 



N} = 



dG 1 



( 7 -l)(2 7 y-l) + (27 2 +7- 2). 



8y 2( 7 -l)(l-2a:) 
The Berwald connection can be obtained as 

dNl 7 



(5.22) 



G n — 



dy 1 — 2x 

Finally, the second KCC invariant or the deviation curvature tensor Pf, defined as 



(5.23) 

reads now 
(5.24) 

pl _ (7 7 -6) 2 x 2 



Pi 



^G 1 

dx 



2G 1 G\ 1 



dNj 
dx 



2(7 - 1)(2 7 2 + 2I7 - 18)a; - 2 7 ( 7 - 1)(2 7 - l)y + 9(7 - l) 2 



4( 7 - l) 2 (l-2x) 2 

Equivalently, this can be written in the variables (u, v) as follows 

x _ (7 7 - 6) V - 4( 7 - 1)(11 7 - 9)u - 27(7 - 1)(2 7 - l)v + 9(7 - l) 2 
(5 ^ 5j n - 4( 7 -l) 2 (l-2 W ) 2 

Evaluating now Pl(u, v) at the steady state (u\, v\), we obtain a rather simple form 
of this function. Equivalently, in the variables (ui, v\) this can be written as follows 

7 2 - 44 7 + 36 



(5.26) 



Pl{ux,v x ) 



4 7 2 



We remark that, for 1 < 7 < 2, we have Pi{u\,v{] < (see Fig. [T2")) . and 
therefore, from the general KCC theory, it follows that the steady point (ui,vi) is 
Jacobi stable without any bifurcation points (see Fig. [T2|) . 

Moreover, we remark that the trajectories of the system (|5.8j) . (|5.9j) are always 
attracted to the steady state (ui,V\) in the interior of the Jacobi stability region. 
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Figure 12. Left: The graph of the function Pl{ui,v{). Right: 
The graph of the function P^(u,v) for 7 = 1.5. The steady state 
(tti,wi) denoted by a point in the graph, is always in the Jacobi 
stability region. 
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Figure 13. The trajectories of the system are always attracted in 
the interior of the Jacobi stability region (we show the graph for 
7 = 1.5). 





In other word, even they start outside the Jacobi stability region, inevitable they 
will end by being Jacobi stable in the final (see Fig. [T3J) • 
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The condition for the Jacobi stability of the star with a linear barotropic equation 
of state, Pl(u,v) < 0, can be written in the initial physical variables m(r),p(r) in 
the form of a restriction imposed on the physical parameters of the stars, and which 
is given by 
(5.27) 



m(r) 



,2 



4(7-l)(ll 7 -l) m(r) 8tt 7 ( 7 - 1) (2 7 - 1) 2 9 ( 7 - l) z ,, 
a 7Z 7^2 P r + — — < °- 



(7 7 -6) z r (7 7 -6) z (7 7 - 

This condition must be satisfied at all points < r < R inside the star. By means 
of some simple algebraical transformations we can rewrite this condition as 
(5.28) 

m(r) 2( 7 -l)(ll 7 -9) | 7 -l / 8tt 7 (2 7 -1) ~ | 4(11 7 ~9) 2 ~ 
r (7 7 - 6) 2 7 7 - 6 y 7 -l PT (7 7 -6) 2 

Due to assuming a linear equation of state for the matter, at the surface of the star 
the matter density and the pressure both vanishes, so that also p(R) = 0. The total 
mass of the star is M = m(R). By evaluating Eq. (|5.28p at the vacuum boundary of 
the star we obtain the following restriction on the mass-radius ratio for barotropic 
stars with a linear equation of state 



(5.29) M < 2 (7-D (117- 9) + T^j j (11 7 - 9) 2 _ 
R (7 7 -6) 2 7 7 - 6 Y (7 7 -6) 2 

When we, for example, choose 7 = 2, then this equations results in 2M/R < 9/8 
which can be combined with our previous assumption M/R > 1/2 (note that this 
condition automatically follow from the quadratic inequality) to give the joint result 

1 2M 9 

(5.30) - < - < -, 

which is a very nice physical result coming out of KCC theory. Note, however, 
that this condition is somewhat weaker than the well-known Buchdahl inequality 
2M/R < 8/9 which implies the absence of horizons in the matter part of the 
spacetime. This cannot be concluded from the KCC approach. 

6. Stability of the vacuum in brane world models 

Brane world models have become a popular modification of Einstein's theory of 
general relativity. In this model one considers a five dimensional spacetime, called 
the bulk spacetime, in which a single four-dimensional brane is embedded. The 
working assumption of these models is that the matter is confined to the brane, 
and only gravity can probe the extra dimension. The action of such a system has 
been formulated by |25[ 126] and a variety of interesting results have been derived 
in the field of brane world gravity. 

When one assumes the four-dimensional manifold to be static and spherically 
symmetric (analogously to the previous Section), then it is possible to derive a set 
of structure equations of the vacuum very similar to Eqs. (|5.1I) - (|5.2[) . The properties 
of the vacuum on the brane are described by the following system of differential 
equations [TTJ US ED H2] 

(6.1) — - — = £ar U, 

dr 
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and 

dU _ {2U + P) [2M + M V - | Ar 3 + a(U + 2P)r 3 ] dP 6P 

\ r t o / 

respectively. In these equations, the quantity Mjj is referred to as the dark mass, 
describing the vacuum gravitational field, exterior to a massive body, in the brane 
world model. U and P are the dark radiation and dark pressure, respectively. A 
is the cosmological constant. The 'matter' source are the projections of the higher 
dimensional gravitational field onto the brane. The brane is apriori in a vacuum 
state. 

To close the system of equations, a supplementary functional relation between 
one of the unknowns U, P and Mjj is needed. Generally, this equation of state is 
given in the form P — P(U). Once this relation is known, Eqs. (|6.1|) - (|6.2[) gi ye a 
full description of the geometrical properties of the vacuum on the brane. 

In the following we will restrict our analysis to the case A = 0. Then the 
system of equations (|6.ip and (|6.2[) can be transformed to an autonomous system 
of differential equations by means of the transformations 

(6.3) u= — + — , v = 3ar 2 U, ir(v) = 3ar 2 P(U(v)), t = lnr. 

r r 

where v and 7r are the 'reduced' dark radiation and pressure, respectively. 

With the use of the new variables given by (|6 . 3[) . Eqs. (|6.ip and (|6.2[) become 

. . du 

(6.4) — = v-u, 

(6.5) * = - (2 - + 7r)[M+( " + 2 - )/3] -2^ + 2^ 27 r. 
dt 1 — q dt 

Eqs. (|6.ip and (|6.2[) , or, equivalently, (|6.4p and (|6.5p . are called the structure 
equations of the vacuum on the brane 11 . In order to close the system of equations 
(I6.4p and (|6.5p an 'equation of state' p = p(fi), relating the reduced dark radiation 
and the dark pressure terms, is needed. 

The structure equations of the vacuum on the brane can be solved exactly in 
two cases, corresponding to some simple equations of state of the dark pressure. 
In the first case we impose the equation of state 2/i + p = 0. From Eq. (|6.5|) 
we immediately obtain p, = Qe~ 29 , while Eq. (|6.4p gives q (9) = —Qe~ 28 + U e~ 9 , 
where Q and Uq = 2GM are arbitrary constants of integration. Therefore we obtain 
the vacuum brane solution 

P = Q_l_ 

2 ~ 3ar 4 

This solution was first obtained in |10j . and therefore corresponds to an equation 
of state of the dark pressure of the form P = —2U. The second case in which the 
vacuum structure equations can be integrated exactly corresponds to the equation 
of state fi + 2p = 0. Then Eq. (|6.5p gives q = 2/3, and the corresponding solution of 
the gravitational field equations on the brane was derived in [11., which corresponds 
to an equation of state of the dark pressure of the form P = —{7/2. 

6.1. Linear stability analysis. Since generally the structure equations of the 
vacuum on the brane cannot be solved exactly, in this Section we shall analyze 
them by using methods from the qualitative analysis of dynamical systems [5] , by 
closely following the approach of [5] . We consider the case in which the dark pressure 



(6-6) U- ,. 
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is proportional to the dark radiation, P = jU, where 7 is an arbitrary constant, 
which can take both positive and negative values. As we have seen, several classes 
of exact solutions of the vacuum gravitational held equations on the brane can be 
described by an equation of state of this form. In the reduced variables u and v the 
linear equation of state is tt = 7U, and the structure equations of the gravitational 
held on the brane have the form 

(6.7) § = 



dv _ 2(1-7) 7 + 2 v [u + 



1+27 



v 



^ 6 ' 8 ) dt (1 + 2 7 ) W l + 2 7 1-u 

Let us firstly analyze the special case where 7 = —1/2. Then, by virtue of 
the second Eq. (|6.8p . we obtain two possible exact solutions, either u = v = 
or it = v — 2/3. The first of these solutions correspond to the vanishing of all 
physical quantities, and therefore we can discard it as unphysical. The second 
case corresponds to an exact solution, which has already been discussed. Let us 
henceforth assume that 7 7^ —1/2. Let us write the dynamical system in the form 

(6.9) § = A£ + B, 



dt 



where we have denoted 



u\ , -1 1 



V ° 2(l- 7 )/(l + 2 7 );' "~^2l2+=r 
The system of equations (|6.9j) has two critical points, Xq = (0,0), and 

3(1-7) 3(1-7) 



(6.11) A- = , — -. — 

v ; 7 V7 2 +7 + 7 7 2 + 7 + 7, 

For 7=1, the two critical points of the system coincide. Depending on the values 
of 7, these points lie in different regions of the phase space plane (u, v). Since the 
term ||-B||/||£|| — > as ||f|| — » 0, the system of equations (|6.9[) can be linearized at 
the critical point Xq. The two eigenvalues of the matrix A are given by n = — 1 and 
r2 = 2(l — 7)/(l + 27), and determine the characteristics of the critical point Xo. 
For 7 £ (—00, — 1/2) U [1, 00) both eigenvalues are negative and unequal. Therefore, 
for such values of 7 the point Xq is an improper asymptotically stable node. 

If 7 € (—1/2, 1), we find one positive and one negative eigenvalue, which corre- 
sponds to an unstable saddle point at the point Xq- Moreover, since the matrix 
dA/d£(Xo) has real non-vanishing eigenvalues, the point Xo is hyperbolic. This im- 
plies that the properties of the linearized system are also valid for the full non-linear 
system near the point Xo. It should be mentioned however, that this first critical 
point is the less interesting one from a physical point of view, since it corresponds 
to the 'trivial' case where both physical variables vanish. 

The structure equations can also be solved exactly for the value 7 = — 2. In 
that case, the non-linear term B in Eq. (|6.9[) identically vanishes, and the system of 
equations becomes a simple linear system of differential equations. For 7 = —2 the 
two eigenvalues of A are given by r\ = — 1 and T2 = —2, and the two corresponding 
eigenvectors are linearly independent. The general solution can be written as follows 

(6.12) £ 7= _ 2 = (uo + wo) (n) e-t + uo ' 



or u 1 
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7 


-oo -0.5 0.67 1 +oo 


r± 


real 


complex 


real 


real 






Re r± < 






r+ 


+ 








r_ 








+ 




Saddle 


Stable spiral 


Stable node 


Saddle 



Table 2. Linear stability of the stable point X~ ( of the static brane 



world vacuum field equations. 



where uq = u(Q) and vq — v(0). One can easily transform this solution back into 
the radial coordinate r form by using t = ln(r), thus obtaining 

(6.13) M 7 = -2 = ^2 , ?7=-2 = — + MO y~ - ^ 

Let us now analyze the qualitative behavior of the second critical point X 1 . To 
do this, one has to Taylor expand the right-hand sides of Eqs. (|6.7[) and (|6.8p around 
Xy and obtain the matrix A which corresponds to the system, linearized around 
X~ r This linearization is again allowed since the resulting non-linear term, n say, 
also satisfies the condition ||n||/||£|| — > as ||£|| — > X T The resulting matrix reads 

(6.14) A = I 3{ _^l 5i _ 4) 

\ (2 +7 )^(l+2 7 ) 

and its two eigenvalues are given by 




,„ ir , -3 - 67 ± v / 167 4 + «7 3 + 132 7 2 - 28 7 - 47 
(6A5) r± = 4 7 2 + 10 7 + 4 • 

For —0.5 < 7 < 0.674865 the argument of the square root becomes negative and 
the eigenvalues complex. Moreover, the values 7 = —1/2 and 7 = — 2 have to be 
excluded, since the eigenvalues in Eq. (|6.14l) are not defined in these cases. However, 
both cases have been treated separately above. 

If 7 € (—00, —1/2) U (l,oo), then X 1 corresponds to an unstable saddle point 
and for 7 e (0.674865, 1) it corresponds to an asymptotically stable improper node. 
More interesting is the parameter range 7 £ (—1/2, 0.674865), where the eigenvalues 
become complex, however, their real parts are negative definite. Hence, for those 
values we find an asymptotically stable spiral point at Xj. Since this point is also 
a hyperbolic point, the described properties are also valid for the non-linear system 
near that point. 

The behavior of the trajectories is shown, for 7 = — 1 and 7 = 0.4, in Fig.Q31 The 
figures show the attracting or repelling character of the steady states, respectively. 
The results of the linear stability analysis of the critical points X 1 are summarized 
in Tabled 



6.2. Jacobi stability analysis. Since from Eq. (|6.7p we can express v as v — 
u + du/dt, upon substitution in Eq. (|6.8p we obtain the following second order 
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q 



Figure 14. Behavior of the trajectories of the structure equations 
of the vacuum on the brane near the critical points for 7 = —1 (left 
figure) and 7 = 0.4 (right figure). 



differential equation 
d 2 u 

+ 



(6.16) dt2 3(l + 27)(l-«) 

+ (4 7 2 + 7 + Yi)u du/dt + (2 7 2 + 57 + 2){du/dtf 



6(7 - l)u + 2( 7 2 + 7 + 7)u 2 + 3(4 7 - l)du 

0. 



which can now be studied by means of KCC theory. 
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By denoting x = u and du/dt = dx/dt = y, Eq. (|6.16|) can be written as 

(6.17) ^ +2 G 1 (x,y) = 0, 
where 

(6.18) G\x,y) = 6 - 1 _ - [6 ( 7 - 1) x + 2 ( 7 2 + 7 + 7) x 2 + 3 (4 7 - 1) y 

+ (4 7 2 + 7 + 13) ay + (2 7 2 + 5 7 + 2) y s 

As a first step in the KCC stability analysis of the vacuum field equations on the 
brane we obtain the nonlinear connection N\ associated to Eq. (|5.19p . and which 
is given by 

, filQ , „i _ dG 1 _ 3 (47 - 1) + (4 7 2 + 7 + 13) + 2 (2 7 2 + 5 7 + 2) y 
( ' 1 ~ 9y ~ 6 (1 + 2 7 ) (1 - x) 

The Berwald connection can be obtained as 

(6 20) G 1 -M- 2 7 2 + 57 + 2 

(6 ' 20) Gll "^T"3(l + 2 7 )(l- a ;)- 

Finally, the second KCC invariant or the deviation curvature tensor Pf, defined 

as 

(6.21) Pi = -2^ - 2Gf 1 Gi 1 + y 9 -§- + 



reads now 



_ 27 - 2 [61 + 57 7 + 4 7 2 (9 + 2 7 )] x + 3 (5 + 7 ) 2 * 2 
(6 - 22) P ^ V) - 12(l-x) 2 (l + 2 7 ) 2 

_ 2(2 + 7 )(l + 2 7 ) (5 + 4 7 )y 
12 (1 - xf (1 + 2 7 ) 2 

Taking into account that x = q and y = /i — g, we obtain P^ in the initial 
variables as 
(6.23) 

27 - 6 (17 + 87 + 2 7 2 ) g + 3 (5 + 7 ) 2 g 2 - 2 (2 + 7) (1 + 2 7 ) (5 + 4 7 ) /x 
1 ^ " 12(l- g ) 2 (l + 2 7 ) 2 • 

Evaluating P* at the critical point X^, given by Eq. (16.111) . we obtain 

(6.24) Pt(X^= ^ + 66 7 -47 

4 (2 + 7 ) 2 (l + 2 7 ) 

The plot of the function Pj 1 (X 7 ) is represented in Fig. [15] 

Using the discussion of Section 16. 1[ our main results on the Linear stability and 
Jacobi stability of the critical point JT 7 of the vacuum field equations in the brane 
world models can be summarized in Table [3] 
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Figure 15. The deviation curvature tensor P\ (X 7 ) as a function 
of 7. 



7 


-00 -C 


.5 0.67 1 +00 


p i( x i) 


+ 




+ 


+ 


Linear stability 


Saddle 


Stable 


Stable 


Saddle 


of X 1 


point 


spiral 


node 


point 


Jacobi stability 


Jacobi 


Jacobi 


Jacobi 


Jacobi 


of X 1 


unstable 


stable 


unstable 


unstable 



Table 3. Linear and Jacobi stability of the stable point X 1 of the 
vacuum field equations on the brane. 



7. Dynamical dark energy models 

Let us consider a dynamically evolving scalar field with exponential potential 
V = Vq cxp(— Xk4>) to model dark energy in a spatially flat Friedmann-Robcrtson- 
Walker (FRW) universe [9]. Moreover, let us assume that the universe also contains 
a dark matter component, modeled as a pressure-less perfect fluid, p 7 = (7 — l)p 7 , 
with 7 = 1. Henceforth we will denote the dark matter density as pdm • The Hubble 
function, describing the dynamics of the expansion of the Universe, is denoted by 
H. 

The evolution equations of this cosmological model are 

(7.1) H 

(7.2) p dm 

(7.3) $ 
subject to the Friedmann constraint 

(7.4) H 2 = ^(p dm + ^ + V). 
where k 2 = 8irG/c 4 . 



j,2\ 



= ^(Pdm + <t> 
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Following [9], we introduce the dimensionless variables u and v defined by 
u := —= — , v := —= — . 

Moreover, let N = ln(a), the evolution equations (|7.ip - (|7.3p can be written as a 
two-dimensional autonomous system of differential equations 

(7.5) ^ = _3 u + A y|, 2 + ^[f + u 2 -, 2 ], 

(7.6) *L = -\Jlu V + l v [i + u 2 -v 2 }, 
while the Friedmann constraint now reads 

2 

K Pdm . 2,2 -i 

where the prime denotes differentiation with respect to N as defined above. 

The critical points of this dynamical system and the results of the linear stability 
analysis are given in Table |H 



Point 


it* 


V* 


Existence 


Stability 


A 








VA 


Saddle point 


B+ 


1 





VA 


Unstable node A < y/6 
Saddle point A > 


B_ 


-f 





VA 


Unstable node A > —y/6 
Saddle point A < — V6 


C 


^3/2 
A 


A 


A 2 >3 


Stable node 3 < A 2 < 24/7 
Stable spiral A 2 > 24/7 


D 


A 

V6 




A 2 <6 


Stable node A 2 < 3 
Saddle point 3 < A 2 < 6 



Table 4. Critical points and linear stability results for the dy- 



namical dark energy models. 



7.1. Lyapunov method. For this particular system under investigation, a suitable 
Lyapunov function is found to be 

(7.7) V(u,v) = (x~x«) 2 + 2(y-y«) 2 

Clearly, V has a global minimum near (u*,v*) and we have 

/ s dV , n ^ T sdu Trt dv . . du . dv 

(7.8 — = (d u V) h (d v V) — = 2 (u- it*) \-Uy-y*) — , 

y ' dN v u 'dN yv ' dN v dN yy y 'dN : 

with du/dN and dv/dN from the autonomous system (|7.5p - (|7.6p . The points A 
and B± are saddle points when analyzed using linear stability analysis and therefore 
we do not expect to find a Lyapunov function for these points. Every function we 
analyzed did not satisfy the required criteria. However, for the points C and D, the 
identified function allows us the establish nonlinear stability. 
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For the point C we analyze the eigenvalues of the Hessian of dV/dN(u, v) at that 
point. The two eigenvalues of Hess(dV/dN(u, v)) are given by 



(7.9) 



9 



-3 - — ± — VA 4 -18A 2 + 90. 
X z X A 



In the region where both eigenvalues are negative, the function dV/dN has a max- 
imum. Since v > 0, this is equivalent with A 2 > 27/8. 
For the point D the two eigenvalues are given by 



(7.10) /i± = -18 + 4A 2 ± ^36 + 6A 2 - A 4 . 

Both eigenvalues are negative if A 2 < 48/17. 

7.2. Jacobi stability analysis. In order to analyze the dynamical system (|7.5|) 
(|7.6j) using the Jacobi method, we first of all have to write this system as one second 
order differential equation of the form 



(7.11) 



x" + 2G\x,y) = 0, 



where x = m, y = v! and the prime denotes differentiation with respect to N. 
Let us start by solving Eq. (|7.5p with respect to v 2 which gives 

(j 12 j 2 _ 2u ' + 3u- 3u 3 

v6A — 3m 

which allows us to compute v' in terms of u and u' . Now, taking the derivative of 
Eq. (|7.5[) and substituting in the expressions of y 2 and yy' gives the second order 
ordinary differential equation F{u" ,u' ,u) = 0. With x = u, y = v! we find that 
G 1 (x,y) is given by 



(7.13) G\x,y) = 



6(3x 2 - 3x 4 + 3xy + y 2 ) 



4(v / 6A - 3u 

+ \V&{-Zx - 3x 3 + 6x 5 - y - 7x 2 y) + \ 2 {6x 2 - 6x 4 + Axy) 

As a first step in the KCC stability analysis of the dark energy model described 
by Eq. (I7.11|) we obtain the nonlinear connection TV 1 = dG 1 /dy which is given by 



(7.14) 



6(3x + 2y) - AV6(1 + 7x 2 ) + 4A 2 



4( V / 6A - 3m) 

The associated Berwald connection G\ x = dN\/dy can be obtained as 

9 



(7.15) 



G n — 



y6A — 3u 

Finally, the second KCC invariant, or the deviation curvature tensor Pf, defined 



as 

(7.16) 

and is found to be 
(7.17) 

Pi(x,y) 



-2 d S-2G^G\ 1+ y d -^ + NlNl 



dx 



8(%/6A-3u) 2 



54a; 2 + Xai(x, y) + \ 2 a 2 (x, y) + \ 3 a3(x, y) + X^a^x, y) 



JACOBI STABILITY . . . GRAVITATION AND COSMOLOGY 



35 



where the four functions at(x,y) are given by 

(7.18) a x {x,y) = 6-y/6(-15a; - 9a; 3 + 12a; 5 - 5y - 7a; 2 ?/), 

(7.19) a 2 (x, y) = 81 + 414a; 2 - 279a; 4 + 168xy, 

(7.20) a 3 (x,y) = 4y/6(-15x + 3a; 3 -2y), 

(7.21) a A (x,y) = 24a; 2 . 

To understand the stability at the critical points A-D, we evaluate P 1 and obtain 



(7.22) 
(7.23) 
(7.24) 
(7.25) 




> 



VA, 



VA ^ 0. 



It is immediately clear that the only point that can be Jacobi stable is in fact 
point C. If A 2 > 24/7 we find that P\{C) < and hence in that case the point is 
stable. Our results for the uncoupled models are summarized in Table [5] 



Point 


A 


Linear 


Jacobi 


Lyapunov 


C 


3 < A 2 < 27/8 
27/8 < A 2 < 24/7 
A 2 > 24/7 


Stable node 
Stable node 
Stable spiral 


unstable 
unstable 
stable 


inconclusive 
asymptotically stable 
asymptotically stable 


D 


A 2 < 48/17 
48/17 < A 2 < 3 
3 < A 2 < 6 


Stable node 
Stable node 
Saddle point 


unstable 
unstable 
unstable 


asymptotically stable 
inconclusive 
inconclusive 


Table 5. Stability analysis of the c 


ynamic dark energy models. 



8. Discussions and final remarks 

In the present paper we reviewed two basic stability analysis methods - (Lya- 
punov's) linear stability analysis and the Jacobi stability analysis, respectively. We 
considered the stability properties of several dynamical systems that play impor- 
tant roles in gravitation and cosmology - Newtonian polytropes, described by the 
Lane-Emden equation, general relativistic fluid spheres, the vacuum gravitational 
field equations of brane world models, and dynamical dark energy models, respec- 
tively. For all cases we used both methods to investigate stability, (Lyapunov) linear 
stability analysis and Jacobi stability analysis, or the KCC theory. The study of 
stability has been done by analyzing the behavior of steady states of the respec- 
tive dynamical system. Linear stability analysis involves the linearization of the 
dynamical system via the Jacobian matrix of a non-linear system, while the KCC 
theory addresses stability of a whole trajectory in a tubular region [23 . 

By using the KCC theory we have been able to derive some physically relevant 
results for each of the considered systems. In the case of the Newtonian polytropes 
from the Jacobi stability condition we obtain a restriction of the ratio of the internal 



3(5 



C. G. BOHMER, T. HARKO, AND S. V. SABAU 



energy of the star and of its gravitational energy in terms of ratio of the density 
and of the mean density of the star. In the case of the general relativistic fluid 
sphere the Jacobi stability condition gives a physical range of the mass-radius ratio 
for general relativistic compact stellar type objects. In the case of the brane world 
models we have shown that the vacuum on the brane is Jacobi unstable for most of 
the values of the parameter 7, with the stability region reduced to a very narrow 
range of 7. For all other values of 7, the vacuum on the brane is unstable, in the 
sense that the trajectories of the structure equations will disperse when approaching 
the origin of the coordinate system. The Jacobi stability analysis also imposes some 
strong restrictions on the parameter A describing the properties of dark energy in 
the dynamic dark energy models. 

It was shown in [23 that the Jacobi stability of a dynamical system can be 
regarded as the robustness of the system to small perturbations of the whole trajec- 
tory. Thus, Jacobi stability is a very conveyable way of regarding the resistance of 
limit cycles to small perturbation of trajectories. On the other hand, we may regard 
the Jacobi stability for other types of dynamical systems (like the ones considered 
in the present paper) as the resistance of a whole trajectory to the onset of chaos 
due to small perturbations of the whole trajectory. This interpretation is based on 
the generally accepted definition of chaos, namely a compact manifold M on which 
the geodesic trajectories deviate exponentially fast. This is obviously related to the 
curvature of the base manifold (see Section [3]). The Jacobi (in) stability is a natural 
generalization of the (in) stability of the geodesic flow on a diffcrcntiablc manifold 
endowed with a metric (Riemannian or Finslerian) to the non-metric setting. In 
other words, we may say that Jacobi unstable trajectories of a dynamical system 
behave chaotically in the sense that after a finite interval of time it would be im- 
possible to distinguish the trajectories that were very near each other at an initial 
moment. 

We have also found that for all the dynamical systems considered in the present 
paper there is a good correlation between the linear stability of the critical points, 
and the robustness of the corresponding trajectory to a small perturbation. Indeed, 
in the case of the gravitational field equations of the vacuum on the brane, for small 
values of the parameter 7 the saddle point is also Jacobi unstable (7 < —0.5 and 
7 > 1), while the stable spiral obtained for —0.5 < 7 < 0.674865 is also robust 
to small perturbations of all trajectories. It is also interesting to remark for the 
same system that for the interval 0.674865 < 7 < 1, the stable node is actually 
Jacobi unstable. In other words, even though the system trajectories are attracted 
by the critical point X 1 one has to be aware of the fact that they are not stable 
to small perturbation of the whole trajectory. This means that one might witness 
chaotic behavior of the system trajectories before they enter a neighborhood of 
Xy. We have here a sort of stability artefact that cannot be found without using 
the powerful method of Jacobi stability analysis. In the present review we have 
presented in detail the main theoretical tools necessary for an in depth analysis of 
the stability properties of dynamical systems that play a fundamental role in our 
understanding of nature. 
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